% Procyclical leverage and crisis probability in a macroeconomic model of bank runs
% Daisuke Ikeda and Hidehiko Matsumoto
% Figure 9

close all
warning off

load('no_policy')
rng('default')

%% Generate random shocks

options = optimset;
T       = 100000;               % simulation periods
t_drop  = 1000;

rng(1)                   % same random shock 
shocks = normrnd(0,1,T+1,1);  % random shock

kv   = zeros(T+1,1);
rbv  = zeros(T+1,1);
nv   = zeros(T+1,1);
xiv  = zeros(T+1,1);
hv   = zeros(T+1,1);
cv   = zeros(T+1,1);
yv   = zeros(T+1,1);
iv   = zeros(T+1,1);
lv   = zeros(T+1,1);
pibv = zeros(T+1,1);
runv = zeros(T+1,1);
rhseev = zeros(T+1,1);
rkhat  = zeros(T+1,1);
rkhats = zeros(T+1,1);
prv  = zeros(T+1,1);
av   = zeros(T+1,1);
erk  = zeros(T+1,1);
x    = zeros(T+1,npol);
x_nonlog    = zeros(T+1,npol);
x_allnonlog = zeros(T+1,npol);
dv   = zeros(T+1,1);
rkhats_next   = zeros(T+1,1);

%% Simulation

zthre = 0;

kv(1)  = stss_k;
rbv(1) = stss_rb;
nv(1)  = stss_n;
lv(1)  = kv(1)/nv(1);
av(1)  = 0;

    for t=2:T+1
        av(t) = rho_a*av(t-1) + sigma_a*shocks(t);
        hv(t)  = ((1-alpha)/psi*exp(av(t))*kv(t-1)^(alpha))^(1/(alpha+1/nu));
        yv(t)  = exp(av(t)) * kv(t-1)^(alpha)*hv(t)^(1-alpha);

        rkhats(t) = log ( rbv(t-1) * (1-1/lv(t-1)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        rkhat(t)  = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
                  + (1+nu)/(alpha*nu+1)*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t-1));
        xiv(t)  = exp(rkhat(t))/exp(rkhats(t));

        K_nonlog  = [log(kv(t-1)) rbv(t-1) log(nv(t-1)) xiv(t)];
        x0_nonlog = (K_nonlog).^expmx;
        x_nonlog(t,:)  = prod(x0_nonlog,2)';
        
        K_allnonlog  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0_allnonlog = (K_allnonlog).^expmx;
        x_allnonlog(t,:)  = prod(x0_allnonlog,2)';

        K  = [kv(t-1) rbv(t-1) nv(t-1) xiv(t)];
        x0 = (log(K)).^expmx;
        x(t,:)  = prod(x0,2)';
        if xiv(t) >= 1
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1)*nv(t-1) - nv(t-1);

            if pibv(t) > 0
                rhseev(t) = exp(x(t,:)*solved_eta1_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta1_rbar);
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0;
            else
                rhseev(t) = exp(x(t,:)*solved_eta3_rhsee);
                rbv(t)    = exp(x(t,:)*solved_eta3_rbar);
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0;
            end

        else
            pibv(t)  = ( exp(rkhat(t)) + 1-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(1+lambda*(1-gamma))*(lv(t-1)-1)*nv(t-1) - nv(t-1);
            nv(t) = n_ss*(1 - nd/100);
            rhseev(t) = exp(x(t,:)*solved_eta2_rhsee);
            rbv(t)    = exp(x(t,:)*solved_eta2_rbar);
            runv(t) = 1;
        end

        cv(t) = rhseev(t)^(-1) + psi*hv(t)^(1+1/nu)/(1+1/nu);
        iv(t) = yv(t) - cv(t);
        kv(t) = (1-delta)*kv(t-1) + iv(t);

        lv(t)  = kv(t)/nv(t);

        erk(t) = log ( alpha*((1-alpha)/psi)^((1-alpha)/(alpha+1/nu)) ) ...
               + (1+nu)/(alpha*nu+1)*rho_a*av(t) - (1-alpha)/(alpha*nu+1)*log(kv(t));

        rbv(t) = fzero('func_Rb',rbv(t),options,lv(t),erk(t),lambda,gamma,delta,sigma_rk);   % deposit interest rate

        rkhats_next(t) = log ( rbv(t) * (1-1/lv(t)) * ( 1+lambda*(1-gamma) ) - (1-delta) );
        z_temp     = (rkhats_next(t) - erk(t)) / sigma_rk;
        cdf = normcdf(z_temp);
        prv(t) = cdf;
        
        dv(t) = kv(t) - nv(t);

    end
    % End of T period simulation

%% pick up runV = 1 after no run some periods

norun_period = 12;

run_num_d  = 0;
duration_q = zeros(1,1);
duration_y = zeros(1,1);
dura_hist  = zeros(1,5);
duration_y_full = zeros(1,1);
dura_hist_full  = zeros(1,25);

for t = t_drop:T-100

    if sum(runv(t-norun_period:t-1)) == 0   % no run in some past periods
        if xiv(t) < 1                       % pick up run
        run_num_d = run_num_d + 1;

            for td = 1:100
                
                % compute credit and output growth for 2 consecutive years
                kvg1 = kv(t+td+4) - kv(t+td  );
                kvg2 = kv(t+td+8) - kv(t+td+4);
                yvg1 = yv(t+td+4) - yv(t+td  );
                yvg2 = yv(t+td+8) - yv(t+td+4);

                if kvg1>0 && kvg2>0 && yvg1>0 && yvg2>0
                    break
                end
                
            end

            duration_q(run_num_d) = td;
            
            for i=1:25
                if duration_q(run_num_d) <= 4*i
                    duration_y_full(run_num_d) = i;
                    dura_hist_full(i) = dura_hist_full(i) + 1;
                    break
                end
            end
            
            if duration_q(run_num_d) <= 4
                duration_y(run_num_d) = 1;
                dura_hist(1) = dura_hist(1) + 1;
            elseif duration_q(run_num_d) <= 8
                duration_y(run_num_d) = 2;
                dura_hist(2) = dura_hist(2) + 1;
            elseif duration_q(run_num_d) <= 12
                duration_y(run_num_d) = 3;
                dura_hist(3) = dura_hist(3) + 1;
            elseif duration_q(run_num_d) <= 16
                duration_y(run_num_d) = 4;
                dura_hist(4) = dura_hist(4) + 1;
            else
                duration_y(run_num_d) = 5;
                dura_hist(5) = dura_hist(5) + 1;
            end
                    
        end
    end
end

close all

mean_dura_q      = mean(duration_q);
mean_dura_y_full = mean(duration_y_full);
dura_sum  = sum(dura_hist);
dura_dist = dura_hist/dura_sum; 

figure('Units','Inches','OuterPosition',[1 0 7 7])
axes1 = axes;
hold(axes1,'on');
hist  = bar(dura_dist,'DisplayName','Model','FaceColor',[0.3 0.75 0.93]);
grid on
set(gca,'FontSize',20);
hold on
x = 1:1:5;
obs = [35, 19, 8, 5, 31]./100;
plot(x,obs,'DisplayName','Data','MarkerFaceColor',[0 0 1],'MarkerSize',10,...
    'Marker','diamond',...
    'LineStyle','none',...
    'Color',[0 0 1]);
xlabel({'Duration (year)'});
title({'Durations of banking crises'});
legend(axes1,'show');

return
